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Abstract. The quantum dissipativc dynamics of a tunneling process through double 
barrier structures is investigated on the basis of a rigorous treatment for the first time. 
We employ a Caldeira-Leggett Hamiltonian with an effective potential calculated self- 
consistently, accounting for the electron distribution. With this Hamiltonian, we use 
the reduced hierarchy equations of motion in the Wigncr space representation to study 
the effects of non-Markovian thermal fiuctuations and dissipation at finite temperature 
in a rigorous manner. Hysteresis, double platcau-like behavior, and self-cxcitcd current 
oscillation are observed in a negative differential resistance (NDR) region of the current- 
voltage curve. We find that while most of the current oscillations decay in time in the 
NDR region, there is a steady oscillation characterized by a tornado-like rotation in 
the Wigner space in the upper plateau of the NDR region. 
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The Caldeira-Leggett (or Brownian) Hamiltonian has been apphed to the 
investigation of quantum dissipative dynamics in several fundamental contexts, 
including quantum tunneling [U [21 13], chemical reactions[l], SQUID rings[5], nonlinear 
optical response[6], and quantum ratchets[7]. Because a complete model of quantum 
dissipative dynamics must treat phenomena that can only be described in real time, 
a great deal of effort has been dedicated to the problem of numerically integrating 
equations of motion derived from the Hamiltonian that describe real-time behavior |8l HI 
[To] . Although such equations are analogous to the classical kinetic equations, which 
have proved to be useful in the study of classical transport phenomena, they are 
difficult to derive in a quantum mechanical framework without approximations and/or 
assumptions. In this paper, we demonstrate that the reduced hierarchy equations of 
motion (HEOM) in the Wigner space representation provide a powerful method to 
study quantum dissipative dynamics in systems subject to non-Markovian and non- 
perturbative thermal fluctuations and dissipation at finite temperature [HI [121 [13]. As 
an example, we employ a model describing the thermal effects in resonant tunneling 
diodes (RTDs) [HI [15] . 

Due to quantum effects, an RTD system exhibits novel negative differential 
resistance (NDR) in its current- voltage (I-V) relation[T6]. Moreover, current 
oscillations [17], plateau-like behavior and hysteresis of the I-V curve have been 
observed in the NDR region[T8]. Theoretically, Frensley found NDR in the I-V curve 
in a numerical computation treating a quantum Liouville equation in the Wigner 
representation that ignored phonon-scattering processes[T9l [20]. Kluksdahl et al. 
incorporated dissipative and self-consistent effects by employing the Poisson-Boltzmann 
equation by adopting a relaxation time approximation and succeeded in modeling the 
experimentally observed hysteresis behavior of the I-V curve [21]. Jensen and Buot 
developed a numerical scheme to treat systems of the same kind and reported that 
current oscillation and plateau-like behavior arise from intrinsic bistability[22[ [23]. 

In the previous studies, dissipative effects of electron transport have been taken 
into account by using the Boltzmann equation. The quantum version of the Boltzmann 
equation, however, employs several approximations and assumptions, and for this 
reason, the validity of that equation in treating quantum dissipative dynamics is 
questionable. We wish to investigate dissipative dynamics in a more quantum 
mechanically rigorous manner. With this aim, we employ a model based on the Caldeira- 
Leggett (or Brownian) Hamiltonian [21 [1], 
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Here, m, p and q are the mass, momentum and position variables of an electron, and 
ruj, pj, Xj and Uj are the mass, momentum, position and frequency variables of the 
jth phonon oscillator mode. The quantities Cj are coefficients that depend on the 
nature of the electron-phonon coupling. Although the validity of the above Hamiltonian 
for describing electron transport phenomena has not yet been investigated thoroughly, 
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because it has a firm quantum mechanical foundation, and because we feel that it is 
necessary to investigate the validity of the previous theoretical results in a fully quantum 
mechanical treatment, we believe that an investigation employing eq.(l) will be helpful 
in constructing a general understanding of electron transport phenomena. 

In eq.(l), U{q,t) is an effective potential for an electron. It can be written 
U{q; t) = UstaticiQ) + UgeifiQ] t), where Ustatici.^) is the static part, determined from the 
conduction band structure, and Useif{,T, t) is the self-consistent part, calculated from 
the electron distribution using the Poisson equation. The heat bath can be characterized 
by the spectral distribution function, defined by J{uj) = J2jcjS{u — Uj)/ {2mjUj), and 
the inverse temperature, /3 = 1/kT, where k is the Boltzmann constant. We assume 
the Drude distribution, given by J{uj) = {m^/n) {•y^u/^'y^ + u"^)), where the constant 
7 represents the width of the spectral distribution of the collective phonon modes 
and is the reciprocal of the correlation time of the noise induced by phonons. The 
parameter ( is related to the electron-phonon coupling strength. We then employ the 
reduced hierarchy equations of motion (HEOM) approach, which can be used to treat 
non-Markovian and non-perturbative system-bath interactions at finite temperature 
with no approximation [T3] . In the HEOM formalism, the reduced density operator 
is expressed in terms of the auxiliary hierarchy density matrix elements, Pj^\..jj(, where 
the indices n and jk arise from the hierarchal expansion of the decay functions e""^* 
and e~'^''^, where z/^ = 2Tik/ (5% is the /cth Matsubara frequency. p!3l Then the 0th 
element is identical to p§})^...fl{t) = Tr^ {ptot(^)}- The HEOM has been used to study 
chemical reactions [I2l [25], linear and nonlinear spectroscopy |26l [27] . exciton and electron 



transfer [281 [291 [301 [3ll [32], and quantum information [331 El]. The HEOM is ideal 
for studying quantum transport systems when we employ the Wigner representation, 
because it allows us to treat continuous systems utilizing open boundary conditions and 
periodic boundary conditions [T2l [20]. 

The equations of motion are written in hierarchical form as follows [T3l [26] : 
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representation is defined by 
t)W{p\q) with U^{p,q) = 
20] . The other operators appearing in eq.(2) 
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If the combinations of n and jk are 



sufficiently large that the condition N = n + Efc=i jfc ^ ^c/ min(7, ui) holds, where ojc 
is the characteristic time scale of electron motion, the infinite hierarchy in eq.(2) can 
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Figure 1. (a) The structure of the RTD. The contact regions are GaAs doped 
2xl0^^cni"^ (42.375 nm), the spacer layers are undoped GaAs (2.825 nm), the barriers 
are undoped AlGaAs (2.825 nm), and the well is undoped GaAs (4.52 nm). (b) The 
averaged potential U{q) for various values of the bias voltage. The potential on the 
emitter side is flat when the biases are 0.270V and 0.322V. Contrastingly, for bias 
voltages in the NDR region, a basin is formed on the emitter side. This basin becomes 
broad for 0.290V. The basins are deepest when the biases are 0.300V and 0.310V, 
corresponding to the lower plateau region of the I-V curve in Fig. 2 (a). The inset 
plots the Fourier components of the current oscillation for 0.290V. 



be truncated at those values of n and jk with neghgible error by setting Wj^J.. j^(t) = 

-{L,^.+^')w^:i^^{tmm- 

We now discuss our numerical simulations using the model described above. As 
the static potential, we employed a double-barrier structure, which models a hetero- 
structure of GaAs sandwiched between two thin AlGaAs layers. The conduction band 
edge consists of a single quantum well bounded by tunneling barriers. The effective 
mass of an electron was assumed to be constant across the device and equal to 0.067mo, 
where mo is the electron mass in vacuum. The size of the position momentum space 
is 100nmx0.704nm~^, which was discretized onto a 356x200 lattice. The structure 
of the RTD is illustrated in Fig. 1 (a). The potential barriers are 0.27eV. We set 
the parameters used in the HEOM as 7 = 12.1THz (7"^ = 8.27fs), C = 72.5GHz 
(C~^ = 13.8ps), and T=300K to create conditions close to those used in previous 
theoretical studies p!9l EH [201 ESI [23]. The depth of the hierarchy and the number of 
Matsubara frequencies were chosen as G {2, 3, 4, 5, 6} and K G {1, 2, 3}, respectively. 
The third-order up- wind and down-wind difference schemes respectively were used in the 
positive and negative p regions to approximate the spatial derivative of the kinetic term, 
— (p/m)d/dq, in order to facilitate implementation of the inflow and outflow boundary 
conditions |1 9 1 [20] . Fourth-order centered difference schemes were employed for the 
other derivatives of p. The fourth-order Runge-Kutta method was used for the time 
evolution. The inflow boundary conditions are set by stipulating Wj'^J..j^{p < 0,q = L) 
and W^'^lj^ip > 0, g = 0) to be given by the equilibrium distribution of a free particle 
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Figure 2. (a) Steady current as a function of applied bias voltage. The black curve 
represents the case in which the bias is increased, and the blue curve represents the 
case in which the bias is decreased, (b) Time evolution of the current for various bias 
voltages in the NDR region. A non-decaying current oscillation is observed in the 
upper plateau region between 0.284V and 0.298V. 

calculated from the HEOM with periodic boundary conditions [12]. Due to fluctuations 
and dissipation, the flow of a wavepacket reaches a steady state even when there exists 
bias voltage. The validity of the boundary conditions was verifled by considering 
several system sizes. At each step, a self-consistent potential Useif{q',t) = —e(f){q]t) 
was calculated using the Poisson equation —e-^(f){q]t) = e[n~^{q) — P(g;t)], where e 
is the dielectric constant [e = 12.85), n^{q) is the doping density {n~^ = 2 x lO^'^cm"^ 
in the doped region), and P{q;t) = f'^oodpWj^] Q{p,q]t) 1271% is the electron density 
calculated from the Wigner distribution. 

We calculated the current-voltage (I-V) characteristics according to the following 
procedure. Under the inflow boundary conditions specifled above, we integrated eq.(2) 
with the effective potential U{q\ t) evaluated iteratively using the Poisson equation 
at zero bias voltage. When the distribution reached the steady state, the current was 
calculated using the distribution \im.i^ooWl^}.Q{p,q]t), and then the flnal distribution 
was used as the initial distribution in the next bias step. We increased the bias from 
O.OOOV to 0.450V and then decreased it to O.OOOV with bias steps of O.OIV and 0.002V in 
the normal and NDR regions, respectively. At each step, we integrated the equations of 
motion until the system reached the steady state distribution. The time periods of the 
integrations were typically 2-4ps in the normal region and 15-40ps in the NDR region. 
The calculated effective potentials (f/(g) = Yimt^oo U{q] t)) for several values of the bias 
voltage in the NDR region are depicted in Fig. 1 (b) . 

The I-V characteristics are presented in Fig. 2 (a). NDR behavior and hysteresis 
in the I-V curve are observed between 0.272V and 0.322V. Moreover, two plateau-like 
structures appear between 0.282V and 0.298V and between 0.300V and 0.318V in the 
case of increasing bias. The time evolution of current in the NDR region for several 
values of the bias voltage is presented in Fig. 2 (b). It is seen that although the 
currents exhibit transient oscillation just after the bias voltage is changed, in all cases 
considered here, other than that of 0.290V, this oscillatory behavior decays after several 
picoseconds. A similar I-V curve was obtained by Jensen and Buot[22]. However, while 
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Figure 3. Time averaged electron density (black solid curve) and effective potential 
(black dashed curve) for (a) 0.290V, (b) 0.300V and (c)0.322V. The red, green, 
blue, orange, and purple curves represent the eigenfunctions in order of increasing 
eigenenergy calculated using the averaged effective potential without the heat bath. 

they observed only one plateau region, our results indicate the existence of two plateau 
regions. In addition, while the oscillation they found does not decay, the oscillation 
found here does besides the cases in the upper plateau. Although multiple plateau 
structures have been observed experimentally. [35| [36| [37] this is the first time that they 
have been reported in a theoretical study. 

In Fig. 1(b), we find a basin-like structure on the emitter side of each potential 
in the NDR region, while no such structure exists in the cases of 0.270V and 0.322V, 
which correspond to the upper and lower boundaries of the NDR regions, as shown in 
Fig. 2(a). These emitter basins behave as a second quantum well and play an important 
role in the appearance of hysteresis, current oscillation and the plateau-like structure of 
the I-V curve [2n [23] . The basin becomes very broad in the case of 0.290V, at which 
value non-decaying current oscillation appears. 

To investigate the current behavior in the NDR region, we plotted the electron 
densities and the averaged effective potentials for several eigenstates. These plots appear 
in Fig. 3: (a) the oscillating upper plateau case (0.290V), (b) a non-oscillating lower 
plateau case (0.300V), and (c) the minimum case (0.322V). There, for the purpose of 
the graph, we employ the time-averaged potential. It should be noted, however, that 
in the simulations, the effective potential varies in time. In the upper plateau case. 
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Figure 4. A snapshot (t=41.4 ps) of the Wigncr distribution function for a bias 
voltage of 0.290V in the case of increasing bias. Current flows into the system from 
the emitter side of the boundary (g^Onm). Then, a part of the current is scattered by 
the emitter side of the barrier. The scattered electron flows in a tornado-like manner 
to central peak in the emitter basin due to dissipation. The shaking motion of the 
effective potential periodically accelerates the component at the central peak to the 
tunneling state, and the current thus exhibits steady oscillation. 

the 3rd (blue) and 4th (orange) eigenstates are the tunnehng states, whereas in the 
lower plateau case the 1st (red) and 2nd (green) eigenstates are the tunneling states. 
Although here we did not observe a steady oscillation, the current oscillation observed 
in the previous studies on the single plateau case corresponds to the lower plateau case, 
since it arises from the ground tunneling state in the basin potential. [38] 

The frequency distribution of the current is displayed in the inset of Fig. 1(b). The 
lower peak there corresponds to transitions between the 5th (purple) and 4th (orange) 
or the 3rd (blue) and 2nd (green) eigenstates, while the higher peak can be attributed 
mainly to transitions between the 3rd (blue) and 1st (red) eigenstates. These results 
indicate that the current oscillation arises from variation of the electron distribution 
in the emitter basin. The upper plateau oscillation is unique in that it arises from 
the excited tunneling states rather than the ground tunneling state. Experimentally 
observed multiple plateau structures [35], |36l [37] may also be explained in terms of this 
origin. 

We display a snapshot of the Wigner distribution for the case of 0.290V in Fig. 
4. The tornado-like distribution centered at g=30nm, appears in the emitter basin 
range between g=10-45nm andp=-0.05— 0.05nm~^. This tornado-like wavepacket rotates 
clockwise, and whenever its tail reaches the maximum value of the momentum, near 
p = O.lnm"^, current flows to the collector side as a result of tunneling. This 
phenomenon is due to the shaking motion of the effective potential, which arises from 
competition between the motion of the electrons and the field that they induce. 

The difference between oscillating and non-oscillating cases in the NDR region can 
be explained on the basis of Fig. 3. In the cases considered in Figs. 3(b) and 3(c), the 
tunneling states are in the bottom of the basin potential. The profiles of the electron 
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distributions in these cases indicate that none of the populated states are the ground 
state, but states close to the emitter voltage. Current oscillation does not occur in 
these cases because the tunneling states are largely unpopulated. If the dissipation were 
larger, the lower tunneling state may be populated. However, due to damping, the 
current oscillation would be suppressed. The lower plateau feature can be explained by 
considering the two adjacent states at the bottom of the basin potential illustrated in 
Fig. 3(b). Because the other states may also result in tunneling, the output voltage 
does not change significantly when the emitter voltage is increased slightly. In the case 
of Fig. 3(c), the basin becomes so narrow that the resonant tunneling state is almost 
completely unpopulated. Once the bias voltage exceeds 0.322V, the basin disappears, 
and the I-V characteristics become normal. 

In summary, we investigated the I-V characteristics and dynamical features of the 
current in the NDR region on the basis of the HEOM approach, employing the Caldeira- 
Leggett Hamiltonian. We find that while most of the current oscillations decay in time 
in the NDR region, there is a steady oscillation characterized by a tornado-like rotation 
in the Wigner space in the upper plateau of the NDR region. To explore the condition 
of current oscillations in the NDR region, calculations for different physical conditions 
are necessary. We will discuss the details of these calculations in a forthcoming paper. 
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